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1  Pressure  -  driven  microchannel  flow 


1.1  Introduction 


As  mentioned  in  our  previous  reports  [1,  2],  Lattice  Boltzmann  (LB)  models  [3,  4,  5,  6] 
are  effective  for  problems  where  both  mesoscopic  dynamics  and  microscopic  statistics 
become  important,  as  in  the  case  of  micro  -  channel  flows  [7,  8,  9,  10,  11,  12,  13].  In 
this  report,  we  further  investigate  the  applications  of  the  two-dimensional  {2D)  thermal 
hnite  difference  Lattice  Boltzmann  (FDLB)  model  introduced  in  [14],  For  convenience, 
a  brief  description  of  this  model  is  provided  below. 

Lattice  Boltzmann  (LB)  models  provide  an  alternative  to  particle-based  models 
like  Molecular  Dynamics  (MD)  or  Direct  Simulation  Monte  Carlo  (DSMC),  as  well  as 
to  numerical  techniques  of  Computational  Fluid  Dynamics  (CFD)  derived  from  the 
principles  of  continuous  media  mechanics.  The  development  of  Lattice  Boltzmann 
models  for  the  investigation  of  flow  phenomena  at  non-negligible  Knudsen  numbers 
[15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26]  has  recently  received  considerable  at¬ 
tention  due  to  the  increasing  interest  for  mastering  gas  flow  at  two  extreme  scales 
[7,  8,  9,  10,  11,  12,  13]:  high-altitude  aerodynamics  and  micro-electro-mechanical  sys¬ 
tems  (MEMS).  Unfortunately,  the  majority  of  LB  models  developed  up  to  date  for  the 
investigation  of  fluid  flow  at  hnite  Knudsen  numbers  still  refer  to  the  isothermal  case. 

In  the  previous  reports  [1,  2]  we  approach  the  non-isothermal  gas  how  in  microchan¬ 
nels  using  the  two-dimensional  hnite  diherence  LB  model  with  multiple  speeds  of  Watari 
and  Tsutahara  [14],  which  allows  the  correct  recovery  of  mass,  momentum  and  energy 
equations  of  a  compressible  huid.  The  model  involves  a  set  of  33  nondimensionalized 
velocities  {eoo,  e^j},  (/c  =  1, . . .  4,  i  =  1, . . .  8) 


Goo  —  0 


^ki  — 


COS 


n{i  —  1) 
4 


sin 


n{i  —  1) 
4 


with  Cfc  G  {1.0,  1.92,  2.99,  4.49}.  The  corresponding  distribution  functions /oo 
fki  =  fki{^,t),  evolve  according  to  the  non-dimensionalized  equations 


(1) 

/oo(x,t). 


dtfki  +  eki-^fki  = -—[fki- fki]  (2) 

T 

Here  we  use  the  density-dependent  relaxation  time  r  =  A/nc,  where  A  is  a  dimensionless 
constant,  n  =  n(x,  t)  is  the  particle  number  density  and  c  is  the  local  average  speed  of 
huid  particles  [27].  The  equilibrium  distribution  functions  in  eqs.  (2) 


fk!  =  fk!i^^t)  =nFkSki 


(3) 


are  expressed  using  the  series  expansion  Ski  =  Ski{9,  u)  up  to  fourth  order  with  respect  to 
the  Cartesian  components  Mq,  (a  =  1,  2)  of  the  huid  velocity  (summation  over  repeated 
Greek  indices  is  understood): 
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The  weight  factors  =  Fk{9)  in  eq.  (3)  depend  on  the  local  temperature  9  =  9{'x,t) 
and  the  speeds  c^,  /c  =  1, . . .  4.  We  refer  the  reader  to  the  literature  [1,  2,  14,  27,  28]  for 
further  details  on  the  thermal  LB  model. 

1.2  Boundary  conditions 

1.2.1  Walls 

The  thermal  LB  model,  as  originally  developed  by  its  authors  [14] ,  uses  the  second  order 
upwind  hnite  difference  scheme  to  solve  the  LB  evolution  equations  (2)  in  the  nodes  of 
a  two-dimensional  square  lattice.  The  boundary  nodes  were  located  on  the  channel 
walls  and  an  extrapolation  method  was  used  to  provide  the  corresponding  values  of  the 
distribution  functions.  This  procedure  forces  the  fluid’s  temperature  and  velocity  in  the 
boundary  nodes  to  equal  the  wall’s  temperature  and  velocity.  Appropriate  boundary 
conditions  [15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26],  many  of  them  based  on  the 
concept  of  diffuse  reflection,  which  dates  back  to  the  time  of  Maxwell  and  Smoluchowski 
[29,  30,  31,  32],  were  already  introduced  in  the  isothermal  LB  models  to  capture  the 
fluid  velocity  slip  near  the  walls  when  the  Knudsen  number  is  no  longer  negligible. 

To  reduce  numerical  errors,  the  MCD  flux  limiter  hnite  difference  scheme  [33,  34, 
35]  is  used  in  our  simulations.  A  hrst  implementation  of  diffuse  rehection  boundary 
conditions  in  the  thermal  LB  model  of  Watari  and  Tsutahara  relies  on  the  redistribution 
(thermalization)  of  the  particle  distribution  functions  in  the  ghost  nodes  outside  the  walls 
[27].  Although  successful  in  capturing  the  temperature  jump  and  slip  velocity  in  Couette 
how,  this  implementation  of  the  dihuse  rehection  boundary  conditions  does  not  allow 
the  huid  density  to  vary  along  the  channel  walls,  as  happens  in  the  case  of  pressure- 
driven  gas  how.  To  overcome  this  drawback,  in  the  previous  report  [2]  we  managed  the 
distribution  functions  to  mix  (redistribute)  themselves  in  wall  nodes  instead  of  ghost 
nodes.  More  precisely,  the  distribution  functions  whose  corresponding  velocities  point 
to  the  wall  in  the  normal  direction  mix  separately  from  those  distribution  functions 
with  corresponding  velocities  orientated  along  the  diagonals  of  the  square  lattice.  For 
convenience,  we  refer  to  the  left  wall  of  the  vertical  channel  in  Figure  1,  where  the  black 
squares  (j  =  1/2,  /)  and  (j  =  1/2, 1  + 1/2),  /  =  0,  1,  . . .,  denote  the  mixing  points  on  the 
wall.  The  corresponding  (unknown)  values  of  the  particle  number  density  are 
and  respectively.  Although  the  mixing  points  may  have  different  temperatures, 

in  this  paper  we  restrict  ourselves  to  the  case  when  the  channel  wall  temperature  is 
constant  {9w‘^’’‘  =  =  9^,  for  all  /). 

In  the  hnite  difference  LB  model  we  need  appropriate  procedures  to  compute  the 
values  of  the  distribution  functions  (/  =  0,  1,  . . .)  dehned  in  the  ghost 

nodes  outside  the  left  wall  of  velocity  u^,  [27].  The  requirement  that  the  distribution 
functions  follow  the  Maxwellian  distribution  law  at  the  mixing  points  on  the  left  wall 
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Figure  1:  Implementation  of  diffuse  reflection  boundary  conditions  in  the  lattice  nodes 
near  the  lower  left  end  of  a  vertical  channel  (the  dashed  line  marks  the  outflow  bound¬ 
ary);  •  -  bulk  nodes,  o  -  boundary  nodes,  □  -  ghost  nodes,  ■  -  wall  nodes  where  the 
distribution  functions  fl\  follow  the  Maxwellian  distribution. 
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reads  (A;  =  1, ...  4): 
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Eqs.  (5),  together  with  the 
wall  in  the  mixing  nodes: 

J2ck 

k 


requirements  that  there  is  no  mass  flux  perpendicular  to  the 


[fli  +  flt\ 


^0,/  ^0,/+! 
Jk2  +  Jk8 
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k 


may  be  solved  to  get  the  values  of  the  distribution  functions  in  the  ghost  nodes  (0,  /) 
and  (0, 1  +  1)  after  each  time  step.  Similar  relations  may  be  used  to  compute  the  values 
of  the  corresponding  distribution  functions  in  the  ghost  nodes  outside  the  right  wall  of 
the  channel. 


1.2.2  Inflow/outflow  boundaries 

The  wall  temperature  and  velocity  are  set  to  6^  =  I  and  =  0,  while  the  pressure  on 
the  inflow  and  outflow  boundaries  are  Pin  and  Pout,  respectively.  At  both  inflow/outflow 
boundaries,  the  x  and  y  components  of  the  fluid  velocity  are  linearly  extrapolated  along 
the  y  direction  to  get  the  corresponding  values  in  the  ghost  nodes.  The  linear  extrapo¬ 
lation  is  used  also  to  define  the  temperature  values  in  the  ghost  nodes  below  the  outflow 
boundary  (Figure  1),  while  the  temperature  in  the  ghost  nodes  above  the  inflow  bound¬ 
ary  is  set  equal  to  the  wall  temperature  6*^,  as  done  by  other  authors  [36,  37,  38].  The 
values  of  the  fluid  density  in  the  inflow/outflow  ghost  nodes  are  calculated  according 
to  the  non-dimensionalized  form  of  the  ideal  gas  equation  of  state  p  =  nO.  Once  the 
values  of  n,  9  and  u  are  known  in  the  inflow/outflow  ghost  nodes,  the  values  of  the  dis¬ 
tribution  functions  in  these  nodes  are  set  to  the  corresponding  values  of  the  equibrium 
distribution  functions,  which  may  be  calculated  according  to  eq.  (3). 

1.3  Lattice  Boltzmann  results:  short  channel 

Figure  2  shows  the  stationary  profiles  of  pressure,  density,  temperature  and  y  component 
of  the  fluid  velocity  u  in  the  stream  direction  down  the  central  line  x  =  L/2  oi  a  short 
channel  with  L  =  0.2,  h  =  0.6,  Pout  =  5  x  10^  and  V  =  Pin/Pout  =  3.  These  profiles 
were  recovered  after  500,000  iterations  with  time  step  St  =  10“^  on  a  square  lattice 
with  21  X  63  nodes  and  spacing  Ss  =  L/21.  We  choose  A  =  10®  to  achieve  the  initial 
value  E^out  =  0.1  of  the  Knudsen  number  Kn  =  A/nL  [27]  on  the  outflow  boundary. 
The  LB  simulation  results  in  Figure  2  may  be  compared  to  DSMC  simulation  results 


particle  number  density  n  pressure 
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Figure  2:  Stationary  profiles  (pressure,  density,  temperature  and  y  component  of  the 
fluid  velocity)  in  the  stream  direction  y  down  the  central  line  a;  =  L/2  of  a  short  channel 
(L  =  0.2,  h  =  3L,  Knout  =  0.1). 
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Figure  2:  (continued)  Stationary  profiles  (pressure,  density,  temperature  and  y  compo¬ 
nent  of  the  fluid  velocity)  in  the  stream  direction  y  down  the  central  line  a;  =  L/2  of  a 
short  channel  {L  =  0.2,  h  =  3L,  Knout  =  0.1). 
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[36,  37,  39,  40],  as  well  as  to  the  BGK-Burnett  solution  of  microchannel  flow  in  the  slip 
regime  [38].  For  this  purpose,  the  prohles  in  Figure  2  should  be  rescaled  taking  into 
account  that  we  used  the  reference  values  po  =  10®,  no  =  10®,  Cg  =  \/2  and  Tq  =  1 
for  pressure,  density,  sound  velocity  and  temperature,  while  the  corresponding  reference 
values  in  [36,  37,  38]  are  Pq  =  6.05  x  10“"^,  no  =  1.21  x  10“®,  Cg  =  0.91  and  Tq  =  1, 
respectively.  A  main  characteristic  of  the  gas  flow  in  the  short  channel,  which  is  well 
captured  by  both  LB  and  DSMC  simulations,  is  the  large  temperature  drop  near  the 
onflow  boundary. 


1.4  Lattice  Boltzmann  results:  long  channel 

Figure  3  shows  the  stationary  profiles  of  temperature  and  pressure  in  the  stream  direc¬ 
tion  down  the  central  line  x  =  Lj^  of  a  long  channel  (21  x  1050  nodes)  with  L  =  0.2, 
h  =  10.  These  prohles  were  recovered  after  2,000,000  iterations,  while  the  other  simu¬ 
lation  paramenters  remained  the  same  as  in  the  case  of  the  short  channel.  Unlike  the 
short  channel  case,  the  temperature  variation  along  the  long  channel  is  less  than  2%. 
The  nonlinear  presure  prohle  along  the  central  line,  due  to  the  competition  between 
the  compressibility  and  rarefection  effects  [7],  is  in  good  agreement  with  the  analytical 
solution  of  Arkilic,  Schmidt  and  Brener  [38,  41]: 

2 

(7) 


p{y) 

Pout 


{QEnout  +  T*)  ~  iP  “  7)(7^  +  1  +  WEiiout) 


y 


The  cross-stream  prohles  in  Figure  4  demonstrate  the  capability  of  the  thermal  LB 
model  to  capture  also  the  non-vanishing  velocity  component  normal  to  the  walls,  due 
to  the  non-uniform  huid  density  in  the  cross-stream  direction  [41].  As  expected,  the  y 
component  of  the  huid  velocity  (whose  cross-stream  prohle  is  parabolic),  as  well  as  the 
slip  velocity  at  the  walls  (Figure  5)  increase  when  approaching  the  outhow  boundary, 
where  pressure  drops  faster  [7,  25,  41]. 

We  investigated  also  the  capability  of  the  thermal  LB  model  with  dihuse  rehection 
boundary  conditions  to  capture  the  minimum  of  the  volumic  how  rate  as  a  function  of 
the  Knudsen  number.  We  follow  [7]  and  express  the  non-dimensionalized  volumic  how 
rate  in  the  middle  of  the  channel  {y  =  L/2)  as 


Q{y)  = 


p  /  n{x,y)uy{x,y)dx 

Jo _ 

-{dp/dy)nL^CR 


n{x,y)uy{x,y)dx 


^  constant  x 


Pin  Pout 


(8) 


where  p  and  n  are  the  mean  huid  pressure  and  density,  dp/dy  ~  (pin  —  Pout)/h  and  we 
took  advantage  of  the  fact  that  the  temperature  is  quite  constant  in  the  channel.  To 
achieve  various  values  of  Kno^,^,  we  changed  the  value  of  Pout  while  preserving  L  =  0.2, 
h  =  10,  P  =  3  and  =  1.  In  Figure  6  we  compare  our  thermal  LB  results  with  those 
obtained  using  the  CLL  model  of  Cercignani,  Lampis  and  Lorenzani  [42],  as  well  as  with 
the  analytical  solutions  in  the  zero-  and  inhnite-Kn  limits,  where  s  =  1.01615  is  the  slip 
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Figure  3:  Stationary  profiles  (temperature  and  pressure)  in  the  stream  direction  y  down 
the  central  line  x  =  L/2  of  a  long  channel  {L  =  0.2,  h  =  SOL,  pin  =  Spout,  W^out  =  0.1). 
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Figure  4:  Stationary  profiles  (temperature,  pressure,  density  and  x  component  of  the 
fluid  velocity)  in  the  cross  -  stream  direction  a;  of  a  long  channel  {L  =  0.2,  h  =  SOL, 
y  Pin  “^Pouti  E-^out  0.1). 
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Figure  4:  (continued)  Stationary  profiles  (temperature,  pressure,  density  and  x  compo¬ 
nent  of  the  fluid  velocity)  in  the  cross  -  stream  direction  a;  of  a  long  channel  {L  =  0.2, 
h  =  SOL,  y  =  L/2,  pin  =  Spout,  E^out  =  0.1). 
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Figure  5:  Stationary  profile  of  the  y  omponent  of  the  fluid  velocity  in  the  long  channel. 


coefficient  [43]: 


Qo  =  l/(6Kn)  +  s  +  (2s^  —  l)Kn  (9a) 

Qoo  ~  log(Kn)  (9b) 

For  Kn  <  1.0,  our  thermal  LB  results  £t  well  to  the  CLL  model  and  the  minimum  of 
the  flow  rate  is  situated  around  Kn  =  0.625.  For  larger  values  of  Kn,  the  thermal  LB 
results  are  situated  between  the  analytical  results  predicted  by  eq.  (9b)  multiplied  by  a 
factor  of  4,  and  those  derived  using  the  CLL  model.  This  behaviour  is  very  similar  to 
the  results  of  Toschi  and  Sued,  recovered  using  an  isothermal  LB  model  for  force-driven 
Poiseuille  flow  [20]. 


mass  flux 
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Knudsen  number  Kn 


Figure  6:  Normalized  mass  flux  Q  as  function  of  Knudsen  number  Kn:  comparison  of 
thermal  LB  simulation  results  with  calculated  values  in  the  CLL  model  [42],  as  well  as 
with  analytical  solutions,  Eqs.(9) 
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2  Force  -  driven  microchannel  flow 

Let  us  consider  the  LB  evolution  equation  with  a  constant  external  force  F.y  =  ma^ 
acting  on  the  fluid  particles  of  mass  m  [20,  23,  33,  34] 

dtfki  +  Gki  ■  V/fci  +  a^{eki^  -  =  -y  [fki  -  fk-]  (10) 

The  boundary  conditions  on  the  left/right  walls  of  the  channel  are  the  same  as  in 
Section  1.2.1  and  periodic  boundary  conditions  are  considered  in  the  vertical  direction. 
The  acceleration  vector  is  orientated  along  the  channel  {a^  =  0,  ay  =  a). 

Simulations  were  done  on  a  500  x  3  lattice  with  spacing  6s  =.  We  hxed  the  non- 
dimensionalized  acceleration  value  a  =  1.0  and  the  time  step  6t  =  1.0  x  10“^.  Two 
values  of  the  mean  fluid  density  were  considered:  n  =  1.0  x  10^  and  n  =  5.0  x  10^. 
The  corresponding  values  of  the  Knudsen  number  are  Kn  =  0.5  and  Kn  =  0.1, 
respectively.  In  the  initial  state  {iter  =  0),  the  fluid  is  at  rest  and  begins  to  flow  under 
the  action  of  the  external  force.  Figures  7  and  8  show  the  evolution  of  the  cross-stream 
density,  temperature,  velocity  and  pressure  prohles  for  both  values  of  n  mentioned  above. 
Although  the  velocity  prohle  is  always  parabolic,  a  striking  feature  is  the  occurence  of 
the  bimodal  temperature  prohle,  as  well  as  the  concave  pressure  prohle.  These  features 
were  observed  also  by  other  authors  [36,  37,  38]  using  Direct  Simulation  Monte  Carlo 
(DSMC)  or  BGK-Burnett  solvers. 

As  done  in  the  case  of  the  pressure  -  driven  how,  in  the  force  -  driven  case  we  checked 
also  the  existence  of  the  minimum  of  the  volumic  how  rate  as  a  function  of  the  Knudsen 
number.  The  normalized  volumic  how  rate  is  shown  in  Figure  9  as  function  of  the 
Knudsen  number.  Unlike  the  pressure  -  driven  case,  the  simulation  results  in  the  force 
-  driven  case  ht  well  to  the  analytical  formula,  Eq.  (9a),  even  at  large  values  of  the 
Knudsen  number. 
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Figure  7:  Evolution  of  density,  temperature,  velocity  and  pressure  profiles  in  force-driven 
microchannel  flow:  n  =  1.0  x  10^,  a  =  1.00,  (5s  =  4.0  x  10“^,  6t  =  1.0  x  10“^. 
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Figure  7:  (continued)  Evolution  of  density,  temperature,  velocity  and  pressure  profiles 
in  force-driven  microchannel  flow:  n  =  1.0  x  10^,  a  =  1.00,  5s  =  4.0  x  10“*^,  5t  = 
1.0  X  IQ-L 
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Figure  8:  Evolution  of  density,  temperature,  velocity  and  pressure  profiles  in  force-driven 
microchannel  flow:  n  =  5.0  x  10^,  a  =  1.00,  6s  =  4.0  x  10“"^,  6t  =  1.0  x  10“^. 
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Figure  8:  (continued)  Evolution  of  density,  temperature,  velocity  and  pressure  profiles 
in  force-driven  microchannel  flow:  n  =  5.0  x  10^,  a  =  1.00,  5s  =  4.0  x  10“*^,  5t  = 
1.0  X  IQ-L 
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Knudsen  number  Kn 


Figure  9:  Normalized  mass  flux  Q  as  function  of  Knudsen  number  Kn:  comparison  of 
force  -  driven  simulation  results  with  calculated  values  in  the  CLL  model  [42],  as  well 
as  with  analytical  solutions,  Eqs.(9) 
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3  Conclusion 

We  investigated  the  two-dimensional  micro  channel  flow  using  a  thermal  Lattice  Boltz¬ 
mann  model  with  appropriate  boundary  conditions.  Two  separate  cases  were  considered: 
the  pressure  -  driven  case  and  the  external  force  -  driven  case.  A  characteristics  of  the 
Lattice  Boltzmann  model  is  the  recovery  of  the  density,  temperature,  velocity  and  pres¬ 
sure  fields  from  the  local  values  of  the  discretized  set  of  distribution  functions,  whose 
evolution  is  governed  by  the  Boltzmann  equation,  which  is  easier  to  manage  than  the 
Navier  -  Stokes  or  Burnett  equations. 

Entrance  and  exit  effects  are  present  in  the  pressure  -  driven  case  and  are  clearly 
seen  especially  in  the  longitudinal  temperature  and  velocity  profiles  when  the  fluid  flows 
in  a  short  channel.  In  long  channels,  the  non-linear  pressure  profile  along  the  center 
line,  the  rarefaction  effects,  as  well  as  the  existence  of  the  so  -  called  Knudsen  minimum 
in  the  plot  of  mass  flow  rate  vs.  Knudsen  number  were  found  to  be  in  good  agreement 
with  literature  results. 
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